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A theory is constructed for dense ionic solutions near charged planar walls that is valid for strong 
inter-ionic correlations. This theory predicts a fluctuation-induced, first-order transition and spon¬ 
taneous charge density ordering at the interface, in the presence of an otherwise disordered bulk 
solution. The surface ordering is driven by applied voltage and results in an anomalous differential 
capacitance, in agreement with recent simulation results and consistent with experimental observa¬ 
tions of a wide array of systems. Explicit forms for the charge density profile and capacitance are 
given. The theory is compared with numerical results for the charge frustrated Ising model, which 
is also found to exhibit a voltage driven first-order transition. 


Recently, experimental observations and molecular 
simulations have suggested a link between long-range 
structural correlations and the electrochemical response 
of a double layer capacitor composed of an ionic liquid 
electrolytepHSj- Specifically, these observations have al¬ 
luded to a possible singular response of the differential ca¬ 
pacitance to changes in the applied electric potential |4]. 
It has been postulated on the basis of molecular simu¬ 
lations that this response results from a competition be¬ 
tween entropic effects of packing and local constraints of 
electric neutrality within the ionic liquid near a planar, 
constant potential electrode[5]. Using general arguments, 
I construct an effective field theory for a symmetric so¬ 
lution of dense ionic media that validates this proposal. 
This theory explains the observed anomalous capacitance 
as a result of a first-order interfacial transition associ¬ 
ated with spontaneous charge ordering at the electrode 
surface. 

The interface between a dense ionic solution and a 
metal electrode has been the subject of much recent 
study, due to the development of ionic liquid-based super¬ 
capacitors that exploit charge separation to create high 
power energy storage devices [6, [ 7 ]. Such concentrated 
electrolyte solutions exhibit inter-ionic correlations that 
render typical mean-field theories developed for dilute 
solutions, such as Gouy-Chapman-Stern theory [8], not 
applicable. Extensions of these theories to account for 
excluded volume have been developed 9] [12], which are 
capable of capturing interfacial layering and a nonmono¬ 
tonic capacitance as a function of applied potential. How¬ 
ever, such extensions typically assume a linearly respond¬ 
ing charge density, which necessitates that their predicted 
response functions are bounded. This contrasts molecu¬ 
lar dynamics simulations of a model of BMIM+PF^ on 
graphite electrodes that indicates a voltage driven struc¬ 
tural transition and divergent capacitance A . Experi¬ 
mental indications of similar emergent long-ranged cor¬ 
relations have been observed in many systems, including 
spontaneous two-dimensional ordering of PF6 on gold [13] 
and free surfaces [14], as well as observations of hystere¬ 
sis upon voltage cycling of CgMIM + Tf 2 N _ on epitaxial 


graphene [3j with observed structural bistability [15]. 

To explain these observations of spontaneous interfa¬ 
cial ordering, I consider the implications of two compet¬ 
ing interactions: 1) short range repulsions that arise from 
packing constraints and can favor spontaneous phase sep¬ 
aration, and 2) long-ranged attractions that arise from 
oppositely charged species and frustrate phase separa¬ 
tion. For a symmetric solution, the lowest-order expan¬ 
sion around a uniform charge density yields an effective 
Hamiltonian, 

^B[</>(r)] = Jr ^ 2 (r) + u(p 4 (r) + y|V0(r)| 2 

where, </>( r) = (p + (r) — p_(r))/p is local excess charge 
density, determined by the relative density of cations, 
p + (r), to anions, p_( r), over the mean liquid density, 
p. While for asymmetric solutions, a cubic order term is 
allowed, for this simplistic case, here it is unallowed by 
symmetry. The effective coulomb coupling Q 2 = (z) 2 /e, 
is screened by the optical contribution to the dielectric, 
e, and z is the magnitude of the charge of the ions. In 
principle, the parameters a, u , and m are dependent on 
temperature and pressure, but here they are taken to all 
be real constants. The parameters u and m are assumed 
to be positive, as is necessary to justify the truncation in 
Eq. [TJ and a is assumed to be negative and small so that 
there is an explicit tendency of the uncharged species to 
demix. This tendency is supported by observations of 
spatial clustering in many simple ionic liquidsjTBl. 

Microscopically, the parameters a,u,m and Q are re¬ 
lated to the screening and bare correlation lengths. This 
relation can be clarified by computing the charge suscep¬ 
tibility in momentum space, 

%(k) = (mk 2 + a + 27rQ 2 /k 2 ) 1 k ^ 0 . (2) 

Within the random phase approximation[IT} [T8] . the De¬ 
bye screening length is identified as 4 = ^/\a\/2nQ 2 , 
and bare correlation length as 4 = y/ m /\ a \ with u = 1, 
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FIG. 1. Geometry and composition of the capacitive cell con¬ 
sidered in this work, namely two ideal parallel plates are sep¬ 
arated in the z direction by a solution of nearly symmetric 
ionic liquid. Under applied potential decaying charge den¬ 
sity waves spontaneously form near the interface. The call 
out contains a snapshot of a typical room temperature ionic 
liquid taken from a molecular dynamics simulation. 


to set the basic energy scale, or equivalently setting 
k^T = 1, where &b is Boltzmann’s constant and T is 
the temperature [19 . For room temperature ionic liquids, 
bare correlation lengths are expected to be on the order 
of the size of the molecule(20l f2llj . For the BMIM+PF^ 
mixture studied in previous simulations, the mean molec¬ 
ular diameter is about 5 A|22|. Typical Debye screening 
lengths in ionic liquids and molten salts are between 1 
and 2 A, due to their low permittivity and large molar 
volume [23] . At the Gaussian level of approximation [24], 
the bulk phase diagram for this model includes a phase 
transition from a disordered phase for 4/4 < 2, where 
the 0(r) = 0 on average, to an ordered, microphase- 
separated state for 4/4 > 2, where 0(r) ^ 0 on aver¬ 
age. This microphase-separated state has a characteristic 
wave vector, q s = l/y/£ s £ c , as can be deduced from the 
maximum in the susceptibility in Eq. [2j and arises gener- 
ically from the competition between interactions acting 
over disparate scales[25]. This specific relationship be¬ 
tween the periodicity of charge density oscillations and 
the correlation and screening lengths has also been ar¬ 
rived at previously [TT]. In that work however, rather 
than postulating an effective Hamiltonian for the charge 
density, Ref. [TT| derived the expression from a modi¬ 
fied Poisson equation. Given typical values of both £ c 
and 4 for ionic liquids are 0(1), at ambient conditions 
such liquids are thus expected to be close to microphase 
separation [26]. It is the proximity of this phase transi¬ 
tion that leads to the anomalous capacitive response, as 


will be shown below. 

To analyze the interfacial behavior of this theory, a 
number of simplifications must be made. First, I con¬ 
sider only the case of an ionic liquid in contact with two 
parallel, identical planar electrodes, a geometry that is 
illustrated in Fig. [l] In this geometry, the system is 
symmetric in the plane parallel to the interface, there¬ 
fore degrees of freedom in the xy plane can be integrated 
out. Second, the separation between the two electrodes is 
assumed to be large compared to 4, so that z can be de¬ 
fined over the domain (0,oo)|27]. The resultant effective 
Hamiltonian per unit area is 

H[<l>( z )] = tfekH- 2 )] - h<j>( 0) + ai</> 2 (0)/2 (3) 

where, Hb, is the Hamiltonian of Eq. [l] evaluated for a 
z-dependent order parameter, (f)(z), divided by the area 
of the system in the xy direction. The parameters, h 
and ai, are phenomenological parameters that account 
for the modulation of the fluctuations at the interface. 
The truncation to second order in <j)(z = 0), restricts the 
analysis to weak interactions between the liquid and the 
elect rode [28]. where a\ and h are both order 1, accommo¬ 
dating non-bond interactions like van der Waals forces, 
and small applied potentials. Previous observations of 
ordering near free interfaces in ionic liquids suggest that 
this approximation is sufficient [14! . 

The field h is related to the chemical potential differ¬ 
ence for positive or negative charge density at the in¬ 
terface. Within the assumption of weak direct surface 
interactions, it is expected to be linearly related to the 
applied potential at the electrode, h oc —\fr[29j. For the 
symmetric system considered in this work, terms of ze¬ 
roth order in T that arise from specific chemical inter¬ 
actions can be neglected. In general, local interactions 
can give rise to terms that shift this dependence by a 
constant. Phenomenologically, neglecting such terms is 
the same as setting the zero of h to the potential of zero 
charge. The parameter a\ describes the ability of the sur¬ 
face to modify the local interactions between anions and 
cations in the electrolyte at the interface and arise due to 
altered packing arrangements near the weakly interacting 
surface. 

While it is not analytically tractable to solve for the 
complete partition function determined by Eq. [3j it can 
be approximated by neglecting fluctuations. The mean- 
field interfacial profile is given by d'H/Scfrfz) = 0 where 
</>(z) is the order parameter profile that minimizes the ef¬ 
fective Hamiltonian. The resultant Euler-Lagrange equa¬ 
tion determines the form of the profile, 

acj)(z) — m\7 2 <j)(z) + 2ttQ 2 J <j>(z')\z — z'\ = 0 . (4) 

For conditions near the bulk phase transition, the solu¬ 
tion away from the boundary is homogeneous and the 
free energy is minimized by cj)(z) = 0. Thus, for small h 
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FIG. 2. Charge density profiles, given by Eq. [ 7 ] for various 
ratios of 4/4 = 1.5 (blue), 0.3 (black) and 0.07 (red) and 
infinitesimal field, h > 0. 


and 4/4 close to 2, it is sufficient to linearize Eq. |4]by 
dropping a term proportional to (j) 3 (z ). 

The intergro-differential equation, together with the 
boundary condition from the surface terms in Eq. [3J 

— h + a\4>{z) — md z (j)(z) = 0 z = 0 , (5) 

the condition that the bulk is homogenous, 


where A = a\/m is the extrapolation length typically 
encountered in surface criticality [31.. Within a self- 
consistent Hartree approximation,[32 ? 33] . the surface or¬ 
der parameter is renormalized, = 0^ F \/r, where T 
is a strictly positive function of a,u and (/sIS,- As the 
external field passes through 0, <$* changes discontinu- 
ously, reflecting the renormalization of the order of the 
phase transition [35]. This discontinuous change of 
signals a first-order interfacial transition, and produces 
long-ranged order in a £ c thick slab parallel to the elec¬ 
trode, commensurate with the amplitude of the charge 
density wave away from the interface remaining finite for 
h —>> 0. This symmetry breaking within the plane of the 
electrode is consistent with the onset of 2d crystallization 
of PF^ ions accompanying the microphase separation ob¬ 
served in molecular simulations [5]. It is also explains 
observations of hysteresis upon electrode charging, [3] as 
nucleating domains of charge oscillations near the elec¬ 
trode surface will require times proportional to 4a c [33] 
that are large at the transition. 

From Poisson’s equation and the charge density in Eq. 
[7] the double layer capacitance can be computed. Specif¬ 
ically, the applied potential is equated to the potential at 
the surface of the electrode by integrating Eq. [ 7 ] twice 
[36 . The capacitance at the potential of zero charge is 
given by a sum of three contributions, 


d z (j){z) = 0 z oo , (6) 

and the constraint of electroneutrality, J r </>(r) = 0, are 
sufficient to determine a unique profile for the charge 
density away from the electrode. The solution of this 
equation with these boundary conditions has the form of 
a damped harmonic function[28]. 

4>{z) = -4 e_z/4 cos ( z Qs + 0) (7) 

cos 6 

where tan# = l/q s £ c , and <4 is the value of the charge 
density at the surface of the electrode. The func¬ 
tional form of Eq. [ 7 ] is routinely used to fit experimen¬ 
tal data[21], and exhibits charge oscillations, or “over¬ 
screening” m, which arise from the finite size of the ions. 
Figure [2] shows representative charge density distribu¬ 
tions for three different values of 4/4- For fixed 4 and 
decreasing 4/4? the profile shows increased layering as 
a consequence of approaching the bulk phase transition. 

While the form of the charge density distribution does 
not change at subsequent levels of approximation, or with 
the incorporation of a cubic term in Eq. [l] the depen¬ 
dence of <fi s on the external field does. Within mean-field 
theory and for symmetric mixtures, the surface order pa¬ 
rameter decreases smoothly as the magnitude of the ex¬ 
ternal field goes to zero, with the functional form 

iMF _ h _ / q\ 

m (2/4 + 1/A) ’ l j 


C/2e = 


1 

4 


l 


+ 4 


dc/) s 

~dh 


h= 0+ 


(9) 


where the first two terms are expected from ideal 
solutions- namely the contribution from standard Gouy- 
Chapman theory that is proportional to the inverse of the 
Debye screening length, and the Stern contribution that 
is proportional to the inverse of the correlation length. 
The last term is proportional to the surface susceptibil¬ 
ity and is singular at an interfacial phase transition as 
anticipated from molecular dynamics results [5 . Equa¬ 
tion 9 establishes the relationship between electrolyte 
fluctuations parallel to the electrode and electrochem¬ 
ical response. Away the surface phase transition, the 
capacitance is found to scale as C oc AT -1 / 2 , as found 
previously [IT]. 

While typical electrostatic calculations anticipate a 
bounded capacitance at finite temperature 11. 37], 
viewed from the perspective of classical statistical me¬ 
chanics the potentially unbounded capacitance is a conse¬ 
quence of diverging correlation lengths encountered at a 
phase transition and their relationship to fluctuation and 
response quantities like the differential capacitance [38]. 
Such collective behavior is reminiscent of charging batter¬ 
ies, where ion intercalation can couple to elastic modes of 
an electrode, resulting in discontinuous changes in accu¬ 
mulated charge[39], For the first-order transition found 
here, the capacitance is expected to diverge at the lo¬ 
cation of the surface ordering transition. For negative 
values of A, the surface can order away from the bulk 
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FIG. 3. Finite size scaling analysis for the lattice model de¬ 
fined in Eq. [To] a) Free energies of the accumulated charge 
density within the layer nearest the electrode for various elec¬ 
trode areas of linear size L — 5 — 16. b,c) Typical config¬ 
urations of the charge density, where grey locates the free 
boundary, si — 0. d) Free energy barrier for inverting the 
surface charge as a function of linear electrode size. Lines are 
guides to the eye and errorbars are the size of the symbols. 


transition [ 30 ] • Specifically, by equating the surface free 
energies of the ordered and disorder phases, coexistence 
conditions are found at h —>> 0 and a critical value of 
A = A*, which is A* = —4(£ c ± 1)/^ in the mean-field 
case and must be solved numerically for the Hartree ap¬ 
proximation, though A* < — (a/ 2 — 1)/4EH. This inter¬ 
facial ordering away from the bulk transition is similar 
to that found in other modulated phases, such as block 
copolymers [41] and ferroelectrics [40, 02]. 

To check the robustness of the approximations em¬ 
ployed above, explicit simulations of a discretized ver¬ 
sion of the system defined by the Hamiltonian in Eq. [T| 
are carried out. Specifically, I study a three-dimensional 
charge frustrated Ising model [43[ defined by 

^[W] = - ± SiS i + |E"( ri i) SiS 3 ( 10 ) 

(fij) 

where q is a reduced charge, {s} denotes the vector of 
Ising-like variables, Si = {0, ±1}, the bracket (i,j) de¬ 
notes a restriction over distinct nearest neighbor pairs, 
Yij is the distance between sites i and j on a three- 
dimensional cubic lattice and v(r) is a Coulomb inter¬ 
action evaluated only on those lattice sites, and asymp¬ 
totically approaches l/|r| as |r| oc. The presence of 
the competing long and short ranged interactions yields 


bulk phase transitions in the same universality class as 
the Hamiltonian in Eq. [l] [44], only with discrete states. 

Consistent with the conditions of typical room tem¬ 
perature ionic liquids, I simulate a confined system in a 
region of the phase diagram where the bulk is disordered, 
but close to the first-order transition into a charge density 
wave phase. The bulk phase diagram for the fully occu¬ 
pied charge neutral case, \si\ = 1, has been determined 
from mean-field theory and explicit simulation [45;. Sim¬ 
ulations are run at q = 1 and k^T = 1.2. The sim¬ 
ulation is embedded in an L 2 x L z volume, with pe¬ 
riodic boundary conditions in the xy plane and nonin¬ 
teracting but conducting boundary conditions along the 
z direction 46 . These boundary conditions explicitly 
generate image charges that have been neglected in the 
theoretical analysis but can screen ion-ion interactions. 
A natural order parameter that distinguishes a sponta¬ 
neously polarized interfacial region from a homogeneous 
bulk is the amount of charge in the layer adjacent to 
the interface, S = * u)> which is extensive in 

the area of the interface, L 2 , and approaches ±L 2 in the 
limit that the interface is filled by positive or negative 
charges. This order parameter is analogous to the ampli¬ 
tude of the charge density wave at the electrode in the 
continuum limit, 0 S . Similarly, a divergent susceptibility 
of S to an external field signals an interfacial phase tran¬ 
sition and is accompanied by an singular capacitance. 

Using an extension of Wang Landau sampling [47]. I 
compute the relative free energy, 

F(S)/k B T = - In (*($-£><<*[* t<])>, (11) 

i 

where (...) denotes ensemble average with fixed ions, 
temperature, and cell volume. The results of these cal¬ 
culations are shown in Fig. [3] Each free energy curve dis¬ 
plays two symmetric minima centered about ±0.7 with 
a large free energy barrier in between. The height of the 
barrier, AF, defined as the difference between the local 
maximum of free energy between —0.5 < S/L 2 < 0.5 
and its global minimum value near S/L 2 = ±0.7 is plot¬ 
ted in Fig. |3jd), for various system sizes. As expected for 
a first-order transition in 2d, this barrier scales linearly 
with L and leads to a voltage dependent capacitance that 
diverges as L 2 . For this symmetric solution, the symme¬ 
try breaking occurs with an infinitesimal applied field. 

Conclusion I have shown how an effective field the¬ 
ory for the coarse-grained properties of ionic liquid 
metal interface can yield insight into the structural 
changes and responses that occur under applied volt¬ 
age. By incorporating additional terms into the order 
parameter expansion [34] or explicit image charges at the 
boundaries 37 , the theory could be generalized to non- 
symmetric mixtures as well as strong direct interactions 
within the interface. By adopting the perspective of the 
capacitance as a fluctuation quantity [38], its relation to 
long-ranged correlations within the ionic liquid becomes 
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transparent. With simple theory and numerical simula¬ 
tion the potential singular behavior of the capacitance at 
a surface phase transition is elucidated. This work high¬ 
lights the importance of explicitly incorporating nonlin¬ 
ear behavior that arises from inter-ionic correlations. 
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